"""
Advanced Thermal Boundary Conditions

Implementation of sophisticated thermal boundary conditions for heat transfer
and multi-physics simulations. Provides support for various boundary condition
types including convective cooling, radiation, conjugate heat transfer, and
time-dependent conditions.

Key features:
- Multiple thermal boundary condition types
- Temperature-dependent boundary conditions
- Conjugate heat transfer interfaces
- Time-dependent thermal loading
- Radiation boundary conditions
- Moving boundary support
- Non-uniform boundary condition distributions

Supported Boundary Condition Types:
1. Dirichlet: Fixed temperature T = T_specified
2. Neumann: Heat flux q" = q_specified  
3. Robin/Convective: h(T_ambient - T_wall) = -k∂T/∂n
4. Radiation: σε(T_ambient⁴ - T_wall⁴) = -k∂T/∂n
5. Mixed: Combination of convection and radiation
6. Conjugate: Coupled solid-fluid heat transfer
7. Periodic: Periodic temperature/flux conditions
8. Time-dependent: Time-varying temperature or heat flux
"""

using LinearAlgebra

# Abstract type for thermal boundary conditions
abstract type ThermalBoundaryCondition end

"""
    DirichletTemperature <: ThermalBoundaryCondition

Fixed temperature boundary condition: T = T_specified

# Fields
- `temperature::Union{Float64,Function}`: Specified temperature (K) or function T(x,y,z,t)
- `time_dependent::Bool`: Whether temperature varies with time
"""
Base.@kwdef struct DirichletTemperature <: ThermalBoundaryCondition
    temperature::Union{Float64,Function}
    time_dependent::Bool = false
end

"""
    NeumannHeatFlux <: ThermalBoundaryCondition

Specified heat flux boundary condition: -k∂T/∂n = q"

# Fields
- `heat_flux::Union{Float64,Function}`: Heat flux (W/m²) or function q(x,y,z,t)
- `time_dependent::Bool`: Whether heat flux varies with time
"""
Base.@kwdef struct NeumannHeatFlux <: ThermalBoundaryCondition
    heat_flux::Union{Float64,Function}        # W/m²
    time_dependent::Bool = false
end

"""
    ConvectiveBoundary <: ThermalBoundaryCondition

Convective heat transfer boundary condition: h(T_ambient - T_wall) = -k∂T/∂n

# Fields
- `heat_transfer_coefficient::Union{Float64,Function}`: Heat transfer coefficient h (W/m²·K)
- `ambient_temperature::Union{Float64,Function}`: Ambient temperature T∞ (K)
- `time_dependent::Bool`: Whether conditions vary with time
"""
Base.@kwdef struct ConvectiveBoundary <: ThermalBoundaryCondition
    heat_transfer_coefficient::Union{Float64,Function}  # W/m²·K
    ambient_temperature::Union{Float64,Function}        # K
    time_dependent::Bool = false
end

"""
    RadiationBoundary <: ThermalBoundaryCondition

Radiation heat transfer boundary condition: σε(T_ambient⁴ - T_wall⁴) = -k∂T/∂n

# Fields
- `emissivity::Float64`: Surface emissivity (0-1)
- `ambient_temperature::Union{Float64,Function}`: Ambient temperature T∞ (K)
- `stefan_boltzmann_constant::Float64`: Stefan-Boltzmann constant (5.67e-8 W/m²·K⁴)
- `time_dependent::Bool`: Whether ambient temperature varies with time
"""
Base.@kwdef struct RadiationBoundary <: ThermalBoundaryCondition
    emissivity::Float64 = 0.8                          # Surface emissivity
    ambient_temperature::Union{Float64,Function}        # K
    stefan_boltzmann_constant::Float64 = 5.67e-8       # W/m²·K⁴
    time_dependent::Bool = false
end

"""
    MixedBoundary <: ThermalBoundaryCondition

Combined convection and radiation boundary condition:
h(T_ambient - T_wall) + σε(T_rad⁴ - T_wall⁴) = -k∂T/∂n

# Fields
- `convection::ConvectiveBoundary`: Convection parameters
- `radiation::RadiationBoundary`: Radiation parameters
- `convection_fraction::Float64`: Weight of convection vs radiation (0-1)
"""
Base.@kwdef struct MixedBoundary <: ThermalBoundaryCondition
    convection::ConvectiveBoundary
    radiation::RadiationBoundary
    convection_fraction::Float64 = 0.7  # 70% convection, 30% radiation
end

"""
    ConjugateInterface <: ThermalBoundaryCondition

Conjugate heat transfer interface between solid and fluid domains.
Ensures temperature and heat flux continuity across the interface.

# Fields
- `solid_conductivity::Float64`: Solid thermal conductivity (W/m·K)
- `fluid_conductivity::Float64`: Fluid thermal conductivity (W/m·K)
- `interface_resistance::Float64`: Thermal contact resistance (m²·K/W)
- `coupled_domain_id::Int`: ID of the coupled domain
"""
Base.@kwdef struct ConjugateInterface <: ThermalBoundaryCondition
    solid_conductivity::Float64
    fluid_conductivity::Float64
    interface_resistance::Float64 = 0.0  # No contact resistance by default
    coupled_domain_id::Int = 1
end

"""
    PeriodicBoundary <: ThermalBoundaryCondition

Periodic boundary condition: T(x₁) = T(x₂) and q"(x₁) = q"(x₂)

# Fields
- `periodic_direction::Symbol`: Direction of periodicity (:x, :y, :z)
- `temperature_offset::Float64`: Temperature difference between periodic boundaries
- `flux_offset::Float64`: Heat flux difference between periodic boundaries
"""
Base.@kwdef struct PeriodicBoundary <: ThermalBoundaryCondition
    periodic_direction::Symbol = :x
    temperature_offset::Float64 = 0.0   # ΔT between periodic boundaries
    flux_offset::Float64 = 0.0          # Δq" between periodic boundaries
end

"""
    MovingBoundary <: ThermalBoundaryCondition

Moving boundary with specified temperature or heat flux.
Useful for moving heat sources, sliding contacts, etc.

# Fields
- `base_condition::ThermalBoundaryCondition`: Base thermal boundary condition
- `velocity::Vector{Float64}`: Boundary velocity vector (m/s)
- `velocity_function::Union{Function,Nothing}`: Time-dependent velocity v(t)
"""
Base.@kwdef struct MovingBoundary <: ThermalBoundaryCondition
    base_condition::ThermalBoundaryCondition
    velocity::Vector{Float64} = [0.0, 0.0, 0.0]
    velocity_function::Union{Function,Nothing} = nothing
end

"""
    ThermalBoundarySystem

System for managing multiple thermal boundary conditions across different
regions of the computational domain.

# Fields
- `boundary_regions::Dict{Symbol,Vector{Int}}`: Boundary node indices for each region
- `boundary_conditions::Dict{Symbol,ThermalBoundaryCondition}`: Boundary conditions for each region
- `thermal_conductivity::Float64`: Material thermal conductivity
- `current_time::Float64`: Current simulation time
- `time_step::Float64`: Current time step size
"""
Base.@kwdef mutable struct ThermalBoundarySystem
    boundary_regions::Dict{Symbol,Vector{Int}} = Dict{Symbol,Vector{Int}}()
    boundary_conditions::Dict{Symbol,ThermalBoundaryCondition} = Dict{Symbol,ThermalBoundaryCondition}()
    thermal_conductivity::Float64 = 0.025  # W/m·K
    current_time::Float64 = 0.0
    time_step::Float64 = 0.001
end

"""
    apply_thermal_boundary_conditions!(T, boundary_system, x, y, z=nothing)

Apply comprehensive thermal boundary conditions to temperature field.

# Arguments
- `T::Array`: Temperature field (modified in-place)
- `boundary_system::ThermalBoundarySystem`: Boundary condition system
- `x::AbstractVector`: x-coordinates
- `y::AbstractVector`: y-coordinates  
- `z::Union{AbstractVector,Nothing}`: z-coordinates (3D only)
"""
function apply_thermal_boundary_conditions!(T::Array, boundary_system::ThermalBoundarySystem,
                                           x::AbstractVector, y::AbstractVector, 
                                           z::Union{AbstractVector,Nothing}=nothing)
    if ndims(T) == 2 && z === nothing
        apply_thermal_boundary_conditions_2d!(T, boundary_system, x, y)
    elseif ndims(T) == 3 && z !== nothing
        apply_thermal_boundary_conditions_3d!(T, boundary_system, x, y, z)
    else
        throw(ArgumentError("Dimension mismatch between temperature field and coordinates"))
    end
end

"""
    apply_thermal_boundary_conditions_2d!(T, boundary_system, x, y)

Apply 2D thermal boundary conditions.
"""
function apply_thermal_boundary_conditions_2d!(T::Array{Float64,2}, boundary_system::ThermalBoundarySystem,
                                              x::AbstractVector, y::AbstractVector)
    nx, ny = size(T)
    dx = length(x) > 1 ? abs(x[2] - x[1]) : 1.0
    dy = length(y) > 1 ? abs(y[2] - y[1]) : 1.0
    
    # Apply boundary conditions for each region
    for (region_name, bc) in boundary_system.boundary_conditions
        if haskey(boundary_system.boundary_regions, region_name)
            apply_boundary_condition_to_region_2d!(T, bc, region_name, 
                                                  boundary_system, x, y, dx, dy)
        end
    end
    
    # Handle standard wall boundaries if not explicitly specified
    if !haskey(boundary_system.boundary_conditions, :left)
        # Default conditions for walls not specified
        apply_default_wall_conditions_2d!(T, boundary_system, x, y, dx, dy)
    end
end

"""
    apply_thermal_boundary_conditions_3d!(T, boundary_system, x, y, z)

Apply 3D thermal boundary conditions.
"""
function apply_thermal_boundary_conditions_3d!(T::Array{Float64,3}, boundary_system::ThermalBoundarySystem,
                                              x::AbstractVector, y::AbstractVector, z::AbstractVector)
    nx, ny, nz = size(T)
    dx = length(x) > 1 ? abs(x[2] - x[1]) : 1.0
    dy = length(y) > 1 ? abs(y[2] - y[1]) : 1.0
    dz = length(z) > 1 ? abs(z[2] - z[1]) : 1.0
    
    # Apply boundary conditions for each region
    for (region_name, bc) in boundary_system.boundary_conditions
        if haskey(boundary_system.boundary_regions, region_name)
            apply_boundary_condition_to_region_3d!(T, bc, region_name, 
                                                  boundary_system, x, y, z, dx, dy, dz)
        end
    end
    
    # Handle default conditions for unspecified boundaries
    apply_default_wall_conditions_3d!(T, boundary_system, x, y, z, dx, dy, dz)
end

"""
    apply_boundary_condition_to_region_2d!(T, bc, region_name, boundary_system, x, y, dx, dy)

Apply specific boundary condition to a 2D region.
"""
function apply_boundary_condition_to_region_2d!(T::Array{Float64,2}, bc::ThermalBoundaryCondition,
                                               region_name::Symbol, boundary_system::ThermalBoundarySystem,
                                               x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64)
    nx, ny = size(T)
    t = boundary_system.current_time
    k = boundary_system.thermal_conductivity
    
    if bc isa DirichletTemperature
        apply_dirichlet_2d!(T, bc, region_name, boundary_system, x, y, t)
        
    elseif bc isa NeumannHeatFlux
        apply_neumann_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)
        
    elseif bc isa ConvectiveBoundary
        apply_convective_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)
        
    elseif bc isa RadiationBoundary
        apply_radiation_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)
        
    elseif bc isa MixedBoundary
        apply_mixed_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)
        
    elseif bc isa ConjugateInterface
        apply_conjugate_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)
        
    elseif bc isa PeriodicBoundary
        apply_periodic_2d!(T, bc, region_name, boundary_system, x, y, dx, dy)
        
    elseif bc isa MovingBoundary
        apply_moving_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)
        
    else
        @warn "Unknown boundary condition type: $(typeof(bc))"
    end
end

"""
    apply_boundary_condition_to_region_3d!(T, bc, region_name, boundary_system, x, y, z, dx, dy, dz)

Apply specific boundary condition to a 3D region.
"""
function apply_boundary_condition_to_region_3d!(T::Array{Float64,3}, bc::ThermalBoundaryCondition,
                                               region_name::Symbol, boundary_system::ThermalBoundarySystem,
                                               x::AbstractVector, y::AbstractVector, z::AbstractVector,
                                               dx::Float64, dy::Float64, dz::Float64)
    # Similar to 2D but extended for 3D - implementation details abbreviated for space
    # In practice, each boundary condition type would have 3D variants
    
    t = boundary_system.current_time
    k = boundary_system.thermal_conductivity
    
    if bc isa DirichletTemperature
        apply_dirichlet_3d!(T, bc, region_name, boundary_system, x, y, z, t)
    elseif bc isa ConvectiveBoundary
        apply_convective_3d!(T, bc, region_name, boundary_system, x, y, z, dx, dy, dz, k, t)
    # ... other 3D boundary condition implementations
    end
end

# Specific boundary condition implementations

"""
    apply_dirichlet_2d!(T, bc, region_name, boundary_system, x, y, t)

Apply Dirichlet (fixed temperature) boundary condition in 2D.
"""
function apply_dirichlet_2d!(T::Array{Float64,2}, bc::DirichletTemperature,
                            region_name::Symbol, boundary_system::ThermalBoundarySystem,
                            x::AbstractVector, y::AbstractVector, t::Float64)
    # Get temperature value (constant or time-dependent)
    if bc.temperature isa Function
        # Temperature is a function T(x,y,t) or T(t)
        temp_func = bc.temperature
    else
        # Constant temperature
        temp_value = bc.temperature
        temp_func = (x, y, t) -> temp_value
    end
    
    # Apply to boundary nodes based on region
    if region_name == :left
        for j in 1:size(T,2)
            T[1,j] = temp_func(x[1], y[j], t)
        end
    elseif region_name == :right
        for j in 1:size(T,2)
            T[end,j] = temp_func(x[end], y[j], t)
        end
    elseif region_name == :bottom
        for i in 1:size(T,1)
            T[i,1] = temp_func(x[i], y[1], t)
        end
    elseif region_name == :top
        for i in 1:size(T,1)
            T[i,end] = temp_func(x[i], y[end], t)
        end
    else
        # Custom region - use boundary_regions if defined
        if haskey(boundary_system.boundary_regions, region_name)
            # Apply to custom boundary nodes (implementation depends on node indexing)
            @warn "Custom region Dirichlet BC not fully implemented for region: $region_name"
        end
    end
end

"""
    apply_dirichlet_3d!(T, bc, region_name, boundary_system, x, y, z, t)

Apply Dirichlet boundary condition in 3D.
"""
function apply_dirichlet_3d!(T::Array{Float64,3}, bc::DirichletTemperature,
                            region_name::Symbol, boundary_system::ThermalBoundarySystem,
                            x::AbstractVector, y::AbstractVector, z::AbstractVector, t::Float64)
    # Get temperature function
    if bc.temperature isa Function
        temp_func = bc.temperature
    else
        temp_value = bc.temperature
        temp_func = (x, y, z, t) -> temp_value
    end
    
    # Apply to 3D boundary faces
    if region_name == :left
        for j in 1:size(T,2), k in 1:size(T,3)
            T[1,j,k] = temp_func(x[1], y[j], z[k], t)
        end
    elseif region_name == :right
        for j in 1:size(T,2), k in 1:size(T,3)
            T[end,j,k] = temp_func(x[end], y[j], z[k], t)
        end
    elseif region_name == :bottom
        for i in 1:size(T,1), k in 1:size(T,3)
            T[i,1,k] = temp_func(x[i], y[1], z[k], t)
        end
    elseif region_name == :top
        for i in 1:size(T,1), k in 1:size(T,3)
            T[i,end,k] = temp_func(x[i], y[end], z[k], t)
        end
    elseif region_name == :front
        for i in 1:size(T,1), j in 1:size(T,2)
            T[i,j,1] = temp_func(x[i], y[j], z[1], t)
        end
    elseif region_name == :back
        for i in 1:size(T,1), j in 1:size(T,2)
            T[i,j,end] = temp_func(x[i], y[j], z[end], t)
        end
    end
end

"""
    apply_neumann_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)

Apply Neumann (heat flux) boundary condition in 2D.
"""
function apply_neumann_2d!(T::Array{Float64,2}, bc::NeumannHeatFlux,
                          region_name::Symbol, boundary_system::ThermalBoundarySystem,
                          x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64,
                          k::Float64, t::Float64)
    # Get heat flux value
    if bc.heat_flux isa Function
        flux_func = bc.heat_flux
    else
        flux_value = bc.heat_flux
        flux_func = (x, y, t) -> flux_value
    end
    
    # Apply Neumann condition: -k∂T/∂n = q"
    if region_name == :left
        for j in 2:(size(T,2)-1)
            q_flux = flux_func(x[1], y[j], t)
            # ∂T/∂x = -q"/k, so T[1,j] = T[2,j] + q"*dx/k
            T[1,j] = T[2,j] + q_flux * dx / k
        end
    elseif region_name == :right
        for j in 2:(size(T,2)-1)
            q_flux = flux_func(x[end], y[j], t)
            # Outward normal is positive x direction
            T[end,j] = T[end-1,j] - q_flux * dx / k
        end
    elseif region_name == :bottom
        for i in 2:(size(T,1)-1)
            q_flux = flux_func(x[i], y[1], t)
            T[i,1] = T[i,2] + q_flux * dy / k
        end
    elseif region_name == :top
        for i in 2:(size(T,1)-1)
            q_flux = flux_func(x[i], y[end], t)
            T[i,end] = T[i,end-1] - q_flux * dy / k
        end
    end
end

"""
    apply_convective_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)

Apply convective boundary condition in 2D.
"""
function apply_convective_2d!(T::Array{Float64,2}, bc::ConvectiveBoundary,
                             region_name::Symbol, boundary_system::ThermalBoundarySystem,
                             x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64,
                             k::Float64, t::Float64)
    # Get convection parameters
    if bc.heat_transfer_coefficient isa Function
        h_func = bc.heat_transfer_coefficient
    else
        h_value = bc.heat_transfer_coefficient
        h_func = (x, y, t) -> h_value
    end
    
    if bc.ambient_temperature isa Function
        T_amb_func = bc.ambient_temperature
    else
        T_amb_value = bc.ambient_temperature
        T_amb_func = (x, y, t) -> T_amb_value
    end
    
    # Apply Robin boundary condition: h(T_amb - T_wall) = -k∂T/∂n
    if region_name == :left
        for j in 2:(size(T,2)-1)
            h = h_func(x[1], y[j], t)
            T_amb = T_amb_func(x[1], y[j], t)
            
            # T_wall = (h*T_amb + (k/dx)*T_interior)/(h + k/dx)
            T[1,j] = (h * T_amb + (k/dx) * T[2,j]) / (h + k/dx)
        end
    elseif region_name == :right
        for j in 2:(size(T,2)-1)
            h = h_func(x[end], y[j], t)
            T_amb = T_amb_func(x[end], y[j], t)
            
            T[end,j] = (h * T_amb + (k/dx) * T[end-1,j]) / (h + k/dx)
        end
    elseif region_name == :bottom
        for i in 2:(size(T,1)-1)
            h = h_func(x[i], y[1], t)
            T_amb = T_amb_func(x[i], y[1], t)
            
            T[i,1] = (h * T_amb + (k/dy) * T[i,2]) / (h + k/dy)
        end
    elseif region_name == :top
        for i in 2:(size(T,1)-1)
            h = h_func(x[i], y[end], t)
            T_amb = T_amb_func(x[i], y[end], t)
            
            T[i,end] = (h * T_amb + (k/dy) * T[i,end-1]) / (h + k/dy)
        end
    end
end

"""
    apply_convective_3d!(T, bc, region_name, boundary_system, x, y, z, dx, dy, dz, k, t)

Apply convective boundary condition in 3D.
"""
function apply_convective_3d!(T::Array{Float64,3}, bc::ConvectiveBoundary,
                             region_name::Symbol, boundary_system::ThermalBoundarySystem,
                             x::AbstractVector, y::AbstractVector, z::AbstractVector,
                             dx::Float64, dy::Float64, dz::Float64, k::Float64, t::Float64)
    # Get convection parameters (similar to 2D)
    if bc.heat_transfer_coefficient isa Function
        h_func = bc.heat_transfer_coefficient
    else
        h_value = bc.heat_transfer_coefficient
        h_func = (x, y, z, t) -> h_value
    end
    
    if bc.ambient_temperature isa Function
        T_amb_func = bc.ambient_temperature
    else
        T_amb_value = bc.ambient_temperature
        T_amb_func = (x, y, z, t) -> T_amb_value
    end
    
    # Apply 3D Robin boundary condition
    if region_name == :left
        for j in 2:(size(T,2)-1), k in 2:(size(T,3)-1)
            h = h_func(x[1], y[j], z[k], t)
            T_amb = T_amb_func(x[1], y[j], z[k], t)
            T[1,j,k] = (h * T_amb + (k/dx) * T[2,j,k]) / (h + k/dx)
        end
    elseif region_name == :right
        for j in 2:(size(T,2)-1), k in 2:(size(T,3)-1)
            h = h_func(x[end], y[j], z[k], t)
            T_amb = T_amb_func(x[end], y[j], z[k], t)
            T[end,j,k] = (h * T_amb + (k/dx) * T[end-1,j,k]) / (h + k/dx)
        end
    # ... similar for other faces (bottom, top, front, back)
    end
end

"""
    apply_radiation_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)

Apply radiation boundary condition in 2D.
"""
function apply_radiation_2d!(T::Array{Float64,2}, bc::RadiationBoundary,
                            region_name::Symbol, boundary_system::ThermalBoundarySystem,
                            x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64,
                            k::Float64, t::Float64)
    # Get ambient temperature
    if bc.ambient_temperature isa Function
        T_amb_func = bc.ambient_temperature
    else
        T_amb_value = bc.ambient_temperature
        T_amb_func = (x, y, t) -> T_amb_value
    end
    
    σ = bc.stefan_boltzmann_constant
    ε = bc.emissivity
    
    # Apply radiation boundary condition: εσ(T_amb⁴ - T_wall⁴) = -k∂T/∂n
    # This is nonlinear and requires iterative solution
    if region_name == :left
        for j in 2:(size(T,2)-1)
            T_amb = T_amb_func(x[1], y[j], t)
            
            # Iterative solution for nonlinear radiation condition
            T_wall = T[1,j]  # Initial guess
            for iter in 1:5  # Newton-Raphson iterations
                # Function: f = εσ(T_amb⁴ - T_wall⁴) + k(T_wall - T_interior)/dx = 0
                T_interior = T[2,j]
                
                f = ε * σ * (T_amb^4 - T_wall^4) + k * (T_wall - T_interior) / dx
                df_dT = -4 * ε * σ * T_wall^3 + k / dx
                
                if abs(df_dT) > 1e-12
                    T_wall_new = T_wall - f / df_dT
                    
                    # Ensure physical temperature range
                    T_wall_new = max(T_wall_new, 50.0)  # Minimum 50K
                    T_wall_new = min(T_wall_new, 2000.0)  # Maximum 2000K
                    
                    if abs(T_wall_new - T_wall) < 1e-6
                        break
                    end
                    T_wall = T_wall_new
                end
            end
            
            T[1,j] = T_wall
        end
    elseif region_name == :right
        # Similar implementation for right boundary
        for j in 2:(size(T,2)-1)
            T_amb = T_amb_func(x[end], y[j], t)
            T_wall = T[end,j]
            
            for iter in 1:5
                T_interior = T[end-1,j]
                f = ε * σ * (T_amb^4 - T_wall^4) - k * (T_wall - T_interior) / dx
                df_dT = -4 * ε * σ * T_wall^3 - k / dx
                
                if abs(df_dT) > 1e-12
                    T_wall_new = T_wall - f / df_dT
                    T_wall_new = max(min(T_wall_new, 2000.0), 50.0)
                    
                    if abs(T_wall_new - T_wall) < 1e-6
                        break
                    end
                    T_wall = T_wall_new
                end
            end
            
            T[end,j] = T_wall
        end
    # ... similar for bottom and top boundaries
    end
end

"""
    apply_mixed_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)

Apply mixed convection-radiation boundary condition in 2D.
"""
function apply_mixed_2d!(T::Array{Float64,2}, bc::MixedBoundary,
                        region_name::Symbol, boundary_system::ThermalBoundarySystem,
                        x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64,
                        k::Float64, t::Float64)
    # Combined boundary condition with weighted contributions
    α = bc.convection_fraction
    
    # Apply convection part
    conv_bc = bc.convection
    T_conv = copy(T)
    apply_convective_2d!(T_conv, conv_bc, region_name, boundary_system, x, y, dx, dy, k, t)
    
    # Apply radiation part
    rad_bc = bc.radiation
    T_rad = copy(T)
    apply_radiation_2d!(T_rad, rad_bc, region_name, boundary_system, x, y, dx, dy, k, t)
    
    # Combine results based on weighting
    if region_name in [:left, :right, :top, :bottom]
        # Get boundary indices based on region
        boundary_indices = get_boundary_indices_2d(region_name, size(T))
        
        for idx in boundary_indices
            i, j = idx
            T[i,j] = α * T_conv[i,j] + (1 - α) * T_rad[i,j]
        end
    end
end

"""
    apply_conjugate_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)

Apply conjugate heat transfer interface condition in 2D.
"""
function apply_conjugate_2d!(T::Array{Float64,2}, bc::ConjugateInterface,
                            region_name::Symbol, boundary_system::ThermalBoundarySystem,
                            x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64,
                            k::Float64, t::Float64)
    # Conjugate heat transfer: ensure temperature and heat flux continuity
    # T_solid = T_fluid and k_solid * ∂T_solid/∂n = k_fluid * ∂T_fluid/∂n
    
    k_solid = bc.solid_conductivity
    k_fluid = bc.fluid_conductivity
    R_contact = bc.interface_resistance
    
    # This is a simplified implementation
    # Full conjugate HT would require solving both solid and fluid domains simultaneously
    
    if region_name == :left  # Interface at left boundary
        for j in 2:(size(T,2)-1)
            # Simplified: assume fluid domain provides boundary condition to solid
            # In practice, would need coupled solution
            
            # With contact resistance: q" = (T_solid - T_fluid)/(R_contact + 1/h_interface)
            # For now, use simple continuity
            if R_contact > 0
                # Account for contact resistance
                T_interface = T[2,j]  # Assume adjacent point represents other domain
                T[1,j] = T_interface / (1 + R_contact * k / dx)
            else
                # Direct continuity
                T[1,j] = T[2,j]
            end
        end
    end
    # ... similar for other boundaries
end

"""
    apply_periodic_2d!(T, bc, region_name, boundary_system, x, y, dx, dy)

Apply periodic boundary condition in 2D.
"""
function apply_periodic_2d!(T::Array{Float64,2}, bc::PeriodicBoundary,
                           region_name::Symbol, boundary_system::ThermalBoundarySystem,
                           x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64)
    # Periodic boundary conditions
    if bc.periodic_direction == :x && region_name in [:left, :right]
        # T(x_left) = T(x_right) + ΔT
        ΔT = bc.temperature_offset
        
        if region_name == :left
            for j in 1:size(T,2)
                T[1,j] = T[end,j] + ΔT
            end
        elseif region_name == :right
            for j in 1:size(T,2)
                T[end,j] = T[1,j] - ΔT
            end
        end
        
    elseif bc.periodic_direction == :y && region_name in [:bottom, :top]
        # T(y_bottom) = T(y_top) + ΔT
        ΔT = bc.temperature_offset
        
        if region_name == :bottom
            for i in 1:size(T,1)
                T[i,1] = T[i,end] + ΔT
            end
        elseif region_name == :top
            for i in 1:size(T,1)
                T[i,end] = T[i,1] - ΔT
            end
        end
    end
end

"""
    apply_moving_2d!(T, bc, region_name, boundary_system, x, y, dx, dy, k, t)

Apply moving boundary condition in 2D.
"""
function apply_moving_2d!(T::Array{Float64,2}, bc::MovingBoundary,
                         region_name::Symbol, boundary_system::ThermalBoundarySystem,
                         x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64,
                         k::Float64, t::Float64)
    # Get velocity (constant or time-dependent)
    if bc.velocity_function !== nothing
        velocity = bc.velocity_function(t)
    else
        velocity = bc.velocity
    end
    
    # Apply base boundary condition first
    base_T = copy(T)
    apply_boundary_condition_to_region_2d!(base_T, bc.base_condition, region_name, 
                                          boundary_system, x, y, dx, dy)
    
    # Modify for moving boundary effects
    # This is a simplified treatment - full moving boundary analysis is complex
    
    # For now, apply the base condition with velocity-dependent modifications
    # In practice, would need to track moving interface and interpolate
    
    if region_name == :left
        for j in 1:size(T,2)
            T[1,j] = base_T[1,j]
            # Could add velocity-dependent corrections here
        end
    end
    # ... similar for other boundaries
end

# Utility functions

"""
    apply_default_wall_conditions_2d!(T, boundary_system, x, y, dx, dy)

Apply default thermal boundary conditions to walls not explicitly specified.
"""
function apply_default_wall_conditions_2d!(T::Array{Float64,2}, boundary_system::ThermalBoundarySystem,
                                          x::AbstractVector, y::AbstractVector, dx::Float64, dy::Float64)
    # Apply adiabatic (zero heat flux) conditions to unspecified boundaries
    default_bc = NeumannHeatFlux(heat_flux=0.0)  # Adiabatic
    
    for wall in [:left, :right, :top, :bottom]
        if !haskey(boundary_system.boundary_conditions, wall)
            # Apply default adiabatic condition
            apply_neumann_2d!(T, default_bc, wall, boundary_system, x, y, dx, dy, 
                            boundary_system.thermal_conductivity, boundary_system.current_time)
        end
    end
end

"""
    apply_default_wall_conditions_3d!(T, boundary_system, x, y, z, dx, dy, dz)

Apply default thermal boundary conditions to 3D walls not explicitly specified.
"""
function apply_default_wall_conditions_3d!(T::Array{Float64,3}, boundary_system::ThermalBoundarySystem,
                                          x::AbstractVector, y::AbstractVector, z::AbstractVector,
                                          dx::Float64, dy::Float64, dz::Float64)
    # Apply adiabatic conditions to unspecified 3D boundaries
    default_bc = NeumannHeatFlux(heat_flux=0.0)
    
    for wall in [:left, :right, :top, :bottom, :front, :back]
        if !haskey(boundary_system.boundary_conditions, wall)
            # Would apply default adiabatic condition for each 3D face
            # Implementation abbreviated for space
        end
    end
end

"""
    get_boundary_indices_2d(region_name, grid_size) -> Vector{Tuple{Int,Int}}

Get grid indices for a boundary region in 2D.
"""
function get_boundary_indices_2d(region_name::Symbol, grid_size::Tuple{Int,Int})
    nx, ny = grid_size
    indices = Tuple{Int,Int}[]
    
    if region_name == :left
        for j in 1:ny
            push!(indices, (1, j))
        end
    elseif region_name == :right
        for j in 1:ny
            push!(indices, (nx, j))
        end
    elseif region_name == :bottom
        for i in 1:nx
            push!(indices, (i, 1))
        end
    elseif region_name == :top
        for i in 1:nx
            push!(indices, (i, ny))
        end
    end
    
    return indices
end

"""
    create_thermal_boundary_system(; thermal_conductivity=0.025) -> ThermalBoundarySystem

Create a new thermal boundary system with default parameters.
"""
function create_thermal_boundary_system(; thermal_conductivity::Float64=0.025)
    return ThermalBoundarySystem(thermal_conductivity=thermal_conductivity)
end

"""
    add_boundary_condition!(boundary_system, region, condition)

Add a thermal boundary condition to a specific region.
"""
function add_boundary_condition!(boundary_system::ThermalBoundarySystem, 
                                region::Symbol, condition::ThermalBoundaryCondition)
    boundary_system.boundary_conditions[region] = condition
end

"""
    update_time!(boundary_system, current_time, time_step)

Update the current time in the boundary system for time-dependent conditions.
"""
function update_time!(boundary_system::ThermalBoundarySystem, current_time::Float64, time_step::Float64)
    boundary_system.current_time = current_time
    boundary_system.time_step = time_step
end